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ABSTRACT 

Shear flows have an important impact on the dynamics in an assortment of different astrophysical objects including accreditation 
discs and stellar interiors. Investigating shear flow instabilities in a polytropic atmosphere provides a fundamental understanding of 
the motion in stellar interiors where turbulent motions, mixing processes, as well as magnetic field generation takes place. Here, a 
linear stability analysis for a fully compressible fluid in a two-dimensional Cartesian geometry is carried out. Our study focuses on 
determining the critical Richardson number for different Mach numbers and the destabilising effects of high thermal diffusion. We 
find that there is a deviation of the predicted stability threshold for moderate Mach number flows along with a significant effect on the 
growth rate of the linear instability for small Peclet numbers. We show that in addition to a Kelvin-Helmholtz instability a Holmboe 
instability can appear and we discuss the implication of this in stellar interiors. 
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1. Introduction 

Understanding the complex dynamic interactions in the interior 
of stars, such as the Sun, is crucially important if we are to de¬ 
velop a physical model of these objects in their entirety. To begin 
to obtain a comprehensive knowledge of the motions in stars it is 
convenient to focus initially on the Sun, which we have the most 
detailed observational evidence for. Helioseismology has shown 
that at the base of the solar convection z one there is a thin region 
of radial she ar called the tachocline (Kos ovichev et al.1119971: 
Tobias 2004). This r egion is belie ved to play a crucial role in the 
solar dynamo (see iSilverslI2008I . and references therein). How¬ 
ever, in spite of the evidence of the existence of the tachocline 
and its importance, there is still a considerable amount of work 
to be undertaken to understand this region using mathematical 
modelling techniques. 

Velocity measurements suggest that the tachocline region is hy- 
drodynamically stable against vertical shear flow dMieschl2007l) . 
However, helioseismology is restricted to large-scale time aver¬ 
aged measurements flChristensen-Dalsgaard & Thotnnsonl2007l) 
and so turbulent motions can be still present on small length and 
time scales. Thus it is very plausible for the tachocline to ap¬ 
pear to be stable, using current helioseismology techniques, but 
actually to be hydrodynamically or magnetohydrodynamically 
unstabl e. Though it is widely assumed that the tachocline is sta¬ 
ble (see iTobias 120041) . ISchatzman et alJ d2000h have shown that 
shear turbulence can appear in a narrow part of the tachocline. 
An unstable tachocline would be significantly different in its dy¬ 
namical interactions from a stable region and so, if we are to un¬ 
derstand the role of this region, for example in the solar dynamo 
we must first understand unstable shear flows in a polytropic at¬ 
mosphere. 

Shear flows occur in a wide variety of natural settings as for ex¬ 
ample in oceanic flows, planetary atmospheres, stars and galactic 
discs. Therefore, there have been a number of previous investi¬ 


gations that examine shear flows in different contexts that can 
help inform our approach to the examination of shear flows in 
stars. 

Previous studies of shear flows have shown that such flows can 
undergo what is known as the Kelvin-Helmholtz (KH) insta¬ 
bility, which develops due to conversion of the available ki¬ 
netic energy of the shear flow into kinetic energy of the dis¬ 
turbances (see lDrazin & Reidll2004 chap. 6). In addition to the 
KH instability, oth er instabilities such as baro clinic instability 
(see lCharnevI 19471) . or the Holmboe instability (iHolmboel 19621) . 
can appear when flows are either rotating or stratified. For our 
study the latter one is of greater interest because it is known that, 
while the KH instability is suppressed by stratification, the more 
slowly growing Holmboe modes becom e dominant with increas¬ 
ing stratification flPeltier & Smythlfl 989l) . 

To study any kind of instability it is convenient to start with in¬ 
vestigating the stability threshold of the system. For the exten¬ 
sively studied KH instability, the necessary criterion for stability 
requires the Richardson number to be greater than 1/4 every¬ 
where in the domain (iMilesl[l961l) . This criterion was derived 
for simplifying assumptions, where the fluid is incompressible, 
inviscid and non-diffusive. However, dropping these simplifica¬ 
tions may alter the stability criterion such that in a system where 
thermal diffusion becomes important, and acts on a smaller time 
scale than bu oyancy , the stability criterion requires a significant 
modification. i Dudisi (119744) and [Zahnj (l 974j) have shown that in 
such systems the product of the Richardson number with the Pe¬ 
clet number is the quantity that indicates stability. The effect of 
thermal diffusion on shear instabilitie s was only studied i n the 
Boussinesq appro ximation by Jones (| 1 9771) . iDudisI (Il974l) and 
more recently by Ligni eres et al. 10 * 9991 ) . such that it is not di¬ 
rectly applicable for stellar interiors where large pressure gradi¬ 
ents have to be considered. In a general fully compressible model 
there is the potential for the stability criterion to be altered as the 
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Mach number is varied. In most stellar regions the Mach number 
is assumed to be small but it can still be potentially significant. 
One example are coronal mass eje ctions where shear flow inst a¬ 
bilities were ob served recently by Ofman & Thompson! (1201 lit . 
iMiczekl (1201 3l) considers a fully compressible fluid in an adia¬ 
batic atmosphere, but the effect of varying all, especially ther¬ 
mal, transport coefficients was not studied. Therefore, this study 
does not capture all relevant effects present in stellar interiors. 
Although, considerable work has been undertaken to examine 
shear flows in a variety of different contexts, no work to exam¬ 
ine shear flows in a polytropic atmosphere has been carried out 
and thus will be what we investigate here. 

In this paper we conduct a linear stability analysis to examine 
both the effect of high thermal diffusion and the effect of com¬ 
pressibility on the onset of shear flow instabilities in a stably 
stratified polytropic atmosphere. While the main focus is on KH 
instabilities the appearance and consequences of a Holmboe like 
instability is investigated. The governing equations are given in 
Sect. [2] along with the numerical method used. Our results are 
presented in Sect0followed by a discussion in Sect. [4] 


where m is the polytropic index. Note, this relation is not valid 
for the perturbed quantities derived in the next section. Due to 
the Schwarzschild criterion the fluid is stable against convection 
if the inequality m > l/(y - 1 ) holds, which is the case for a 
polytropic index of m >1.5. Only stable stratified atmospheres 
will be considered throughout the paper. 

Boundary conditions, at the top and bottom of the of the domain, 
are impermeable and stress-free velocity i.e. 

dll, du v 

u- = —— = —- =0 at z — 0 and z = 1 (5) 

dz dz 

and fixed temperature at the top and bottom: 

T = 1 at z = 0 and T = 1 + 0 at z = 1 ■ ( 6 ) 

To include thermal effects it is necessary to choose the back¬ 
ground temperature in such a way, that it is a stationary solution 
of the heat equation or remains quasi-stationary on time scales 
larger than the thermal diffusion time scale. This results in a tem¬ 
perature and density profile of the form: 


2. Model 


2.1. Governing equations, boundary conditions and 
background state. 


We consider a compressible fluid in a Cartesian domain bounded 
at z = 0 and z = 1 and periodic in x and y directions. The fluid 
is assumed to be an ideal gas with constant dynamic viscosity, 
fj, constant thermal conductivity, k, constant heat capacities c p at 
constant pressure and c v at constant volume. The equations we 
consider, in non-dimensional form, are 


dp 

dt 

d(pu) 

dt 

dT 

~dt 


-V- (pu) 


<rC k IV u + -V(V-u)j - V- (puu) 
- V /7 + 9(m + l)pz 

C kw(y ~ !) | t| 2 + rQ v 2 T 

2P P 

-V- (T u) - (y-2)TV-u 


( 1 ) 

( 2 ) 

(3) 


where p is the density, u the velocity field, T the temperature, 
6 denotes the temperature gradient, and p is the pressure. The 
dimensionless Prandtl number, cr = pc p /K , is the ratio of viscos¬ 
ity to thermal conductivity, the thermal dissipation parameter is 
defined as C k = Kt/(poc p d 2 ) and y = c p /c\ denotes the adiabatic 
index. The strain rate tensor has the form 


T 


ij 


du j duj 2 du k 
dxi dxj '' 3 dx k 


In the dimensionless equations above, all lengths have been 
scaled with the domain’s depth d. Recasting the temperature and 
density in units of T, and p r , the temperature and density at the 
top of the layer, and taking the sound-crossing time, which is 
given as t — d/[(c p - c v )T t ] 1 ' 2 , as the fundamental time it fol¬ 
lows that the pressure p is given in units of p, = ( c p - c v )p t T, 
and the velocity field is given in units of acoustic wave velocity. 
For the background state we assume a polytropic relation be¬ 
tween pressure and density such that the pressure is a function 
of density only i.e. 

Pip) p^ 1+ ") (4) 


T(z) = T,{\ + 6z) (7) 

p(z) = p, (1 + 0z) m ( 8 ) 

where 6 is the dimensionless temperature difference between the 
upper and lower boundaries of the domain. These equations form 
only an equilibrium state if the fluid is at rest or viscous heating 
is negligible. The background velocity profile takes the form 


u(z) = Uq tanh ^ 5 j (9) 

with a shear amplitude Uq and a scaling factor 1 /L„ that con¬ 
trols the width of the shear profile. The boundary conditions 
introduced in equation ( 0 ) and © restrict the shear profile to 
values of L„ which will result in a small enough value of the z- 
derivative at the boundaries. For the static state the temperature 
and density profiles are taken as in equation 0 and equation 
©, respectively. Then, the equilibrium state is characterised by 
Uo(z) = O(z), 0,0) r , 7o(z), poiz) andp 0 (z). Selected background 
profiles for temperature, density and velocity are shown in Fig. 

m 

2.2. Formulation of the Eigenvalue problem 

In a diffusive model there are a number of different time scales 
including the time scale associated with the shear dynamics 
4 = L u /(Uq), the time scale for buoyancy tg = 1/N(z), where 
N(z.) 2 is the Brunt-Vaisala frequency, and the time scale for ther¬ 
mal diffusion 4 = L 2 / Ck . In this paper, we focus on the regime 
where the viscous time scale, 4 = L 2 j p, is much greater than 
any other time scales. This allows us to neglect viscous heat¬ 
ing, which corresponds to the first term on the right hand side in 
equation (0. In addition, the shear flow given by equation 0 is 
in equilibrium only if 4 is much greater than the instability time 
scale, which we verify a posteriori. When 4 becomes compara¬ 
ble with other relevant time scales, the background shear flow is 
not in equilibrium and our analysis would be inappropriate for 
this case. 

We perturb each quantity that appears in equations 0 - 0 such 
that / = /o + 6f and 

6f{x,y,z,t) = f(z)exp(ikx + ily + (t), (10) 
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Fig. 1. Plots of the typical background density, temperature and shear 
flow profiles, a) Temperature profiles for different 9, that were used, b) 
Density profiles for 9 = 2 but three different polytropic indices m. For 
all indices m the atmosphere is stably stratified, c) Shear flow profiles 
with the smallest and largest characteristic length L„ used are shown. 


where 1 e 1 and l e R. are the horizontal wave numbers, and 
C, = £ r + i£j e C, where gives the growth rate of the linear 
instability. 

Note, the equations for the perturbed quantities that we obtain 
do not inherit the same symmetry properties as the well known 
Taylor-Goldstein equation (e.g.. lMilesi(il961l) ). where taking the 
complex conjugate of the eigenfunction and eigenvalue leads to 
the same equation. This symmetry is broken in our set of equa¬ 
tions, because there are still terms linear in k and /. Therefore, 
for our eigenvalue problem there do not necessarily exist two 
complex conjugated solutions where one is decaying and one is 
a growing solution. 

For our initial set of equations a Squire transformation exists 
to transform the three-dimensional problem to a corresponding 
two-dimensional one. The transformation can be written as: 


k 2 = 

k 2 + l 2 

kSu = kdu + ISv 5w — 5w 


U 0 - 

U 0 

k 6 (in + 1) 

c - c op — — - - Sp 

k Q(fh + 1) 


dT = 

k 

-ST 

k 

- i 2 

To - ¥ _To 


po = 

kC k 

WT" 

£ 

o 1 o 

a 

II 

°o 

(11) 


set of equations 


&P 

(Pqu 


£p 0 w 


£6T 


-ikupo - —(wpo) - ikdpUo 
dz 


-ik (po6T + T 0 6p) -p 0 w 


dUo 

dz 


ikpuuUo 


r i -iA d ~ u 1 -u dw 


(12) 


(13) 


dz 


( po6T + Todp) - ikUopow + 6(m + l)6p 


4 d 2 w Id 




C k y / d 2 5T 
Po \ dz 2 


k 2 5T - ikU 0 ST 


dT 0 


dw 


- (y - 1) Tq I iku + — - w——, 


dz 


dz 


(14) 


(15) 


where a similar set of eq uations was derived for a problem in¬ 
cluding magnetic fields bv lTobias & Hughesl(l2004i) . Our system 
is characterised by six parameters m , (9, <x, C/ ( , Uq and L u . Equa¬ 
tions ® - ( IT5] > are numerically solved on a one dimensional 
grid in z-direction that is discretised uniformly, this method is 
adapted from the method used bv lFavier et al.l (120121) . Recasting 
the set of differential equations into the form 


ff = Af, 


(16) 


where the matrix A contains the finite difference coefficients ap¬ 
plied to the discretised eigenfunctions f = ( 5p , u, v, w, 6 T) t , re¬ 
duces the problem to a matrix equation. To find the eigenval¬ 
ues and vectors the Schur factorisation is used dAnderson et al.l 
1 19991 ). For the computation of the relevant coefficients in A, 
a central fourth-order finite differences scheme was used. We 
search for the eigenvector solutions with the greatest real part 
of the eigenvalue £ and where the vertical velocity eigenvector, 
w, vanishes at the boundaries. Ultimately, we aim at undertaking 
non-linear simulations with a pseudo-spectral code and viscosity 
will be mandatory in that case. Therefore, most of the computa¬ 
tions will consider a viscous fluid. 


3. Results 

In this section we focus on a number of key areas of interest. 
First, we present the change of the stability threshold while the 
Mach number is varied, which correspond to a continuous tran¬ 
sition between an incompressible and a compressible fluid. The 
effect of compressibility is separately investigated in a weakly 
thermally stratified and a strongly thermally stratified atmo¬ 
sphere Sect. 13.11 Later in Sect. 13.21 the growth rates of the lin¬ 
ear shear instability together with the critical Ri for different Pe- 
clet numbers are compared and the effect on the stability against 
buoyancy is discussed. In Sect. l3.2.H the effect of different poly¬ 
tropic indices on the instability is addressed and the possibility 
of a Holmboe like instability is investigated in Sect. 13.31 

3. 1. The effect of varying the Mach number on the instability 
threshold 


Having also checked numerically that indeed for a certain wave 
number k the growth rate £ r decreases with increasing /, we can 
set l — 0 without loss of generality for our following computa¬ 
tions. Denoting bu = (u, v, w), we obtain this linearised coupled 


As the Richardson criterion is based on simple energetic argu¬ 
ments and does not take compressibility into account, clarifica¬ 
tion is needed to determine whether compressibility affects the 
stability of a shear flow. Therefore, in this section we focus on 
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the stability threshold for different Mach numbers in a viscous 
(we consider cr — 1.0 and C* = 10 6 ) and stably stratified fluid 
with m = 1.6. Although stellar interiors have typically low Mach 
numbers, especially at the base of the convection zone, gener¬ 
ally moderate to high Mach numbers can appear at the surface 
and in other astrophysical objects. Thus, investigating the conse¬ 
quences of moderate Mach numbers on a shear flow is of general 
interest. In the following we refer to the Mach number, M(z), as 

M(z) = ~^=- (17) 

VTTft 


This is the consequence of our previous definition, where veloc¬ 
ity is given in units of the sound speed that is computed at the 
top of our domain. As the inflexion point of our shear flow is at 
Z = 0.5, and the sound speed varies with temperature, it is nec¬ 
essary to com pute the actual M ach number, M, at z = 0.5 , _ 

According to ISchochetl (1 19941) and iGuillard & Murronel (120041) 
the solutions of the compressible Euler equations reduce to the 
solutions of the incompressible Euler equations in the low Mach 
number limit. Thus, varying the Mach number allows to investi¬ 
gate the validity of the Richardson criterion for low to moderate 
Mach numbers (0.02 < M < 0.15). 

We make use of the general definition of the Brunt-Vaisala fre¬ 
quency given by 




(18) 


where T = (P t /P) 1 1/7 is the potential temperature, to define the 
local Richardson number as 


Pin 


mm 

0<z<l 


N(z} 


7 du(z) \ 2 ' 


= min 

0<z<l 


l dz ) , 

e 2 (m + l)(^i - III) 

( l + fe )( t / o _^)) 2 


(19) 


where the derivative of the background velocity profile with re¬ 
spect to z corresponds to a local turnover rate of the shear. 

To find the critical Richardson number, Ri c , we solve the eigen¬ 
value problem for a small Ri, while varying the wave number, 
k, between 0 and 1 /L u to find the most unstable mode k max . For 
large, but finite, Reynolds number the system is assumed to be 
stable if the growth rate, /V, is zero for all wave numbers or the 
time scale for the instability, t( = l/£ r , compared to the viscous 
time scale, t fI , is of the same order. 

A detailed survey for 8 = 2, a weakly stratified atmosphere, 
and for 8 = 10, a strongly stratified atmosphere, reveals that the 
critical Richardson number decreases for Mach numbers greater 
than 0.08, which can be seen in Fig. [2] For the majority of in¬ 
compressible, and weakly compressible Mach numbers the crit¬ 
ical Richardson number does not significantly deviate from the 
well known 1/4 threshold for stability. In the case of a weakly 
stratified atmosphere the critical Richardson number decreases 
rapidly below 0.1 for M « 0.14. In a strongly stratified fluid 
we find qualitatively the same behaviour as for the weakly strati¬ 
fied case but the critical Richardson number does not drop below 
0.2 for M « 0.14. The shift of the stability threshold, in both a 
strong and a weakly stratified fluid, towards smaller Richard¬ 
son numbers for moderate Mach numbers indicates a stabilis¬ 
ing effect of compressibility on the KH instability. Though not 
only focusing on KH like instabilities, previous investigations of 
compressible shear flows with uniform temperature, and without 
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Fig. 2. In both plots the critical Richardson number is found in a viscous 
fluid for different Mach numbers, M. The corresponding wave numbers 
k max of the most unstable mode at the onset of instability are plotted, the 
wave number is normalised by the inverse of the characteristic length 
1/L„. The horizontal line in both plots correspond to Ri = 0.25. In a) 
the atmosphere is weakly stratified with 0 = 2 whereas in b) 8 = 10 
which corresponds to a strongly stratified atmosphere. 


gravit y, for high Mach num bers revealed similar result s. For ex¬ 
ample iBIumen et al.l (1 1 9751) and lDrazin & Davevl (Il977l) showed 
that for M —> V2 the system becomes stable. In this system two 
types of unstable modes are present for Mach numbers greater 
than 0.94, which are stationary modes and travelling modes. 

In the case of KH like instabilities one explanation for the sta¬ 
bilising effect for Mach numbers greater than 0.08 is as fol¬ 
lows. The Richardson criterion uses simple energetic arguments, 
where two neighbouring fluid parcels are exchanged. The den¬ 
sity of these parcels remains constant for an incompressible 
fluid such that only the change in velocity (at different heights) 
changes the kinetic energy, A E^ n , and changes in the potential 
energy, A E pot , are solely due to the changes in position of the 
fluid parcel. While for a compressible fluid the density will de¬ 
crease to a certain amount when a fluid parcel is moved up adia- 
batically, such that a part of the kinetic energy is converted by the 
process of expansion. Therefore, to reach the instability thresh¬ 
old more kinetic energy is needed which requires a greater ve¬ 
locity gradient. 

Different temperature gradients have a non-trivial effect on 
the effective Mach number throughout our domain, where M 
changes according to equation (flTt . In the limit of weak thermal 
stratification M remains almost constant whereas in a strongly 
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Fig. 3. Eigenfunctions for a fixed M = 0.114 and a fixed characteristic 
length scale L u = 0.09 are found for two different 9, where 9 = 2 for 
the plots d, e, f and 9 = 10 for top plots a, b, c. In both cases the 
eigenfunctions for the most unstable mode are shown. 



stratified atmosphere it changes significantly and generates an 
asymmetry between the regions above and below the shear flow. 
The observed asymmetry changes the eigenfunctions found for 
a fixed M{z = 0.5) and L u , which are shown in Fig. [3] Note 
that because M(z = 0.5) and L u are fixed, increasing 9 is equiv¬ 
alent to increasing the Richardson number (see equation (ITOl o. 
The asymmetry with respect to the mid-plane is clearly visible 
in the temperature and vertical velocity eigenfunctions for the 
strongly-stratified case 6 — 10 (see Fig. Q]a and b) whereas the 
same eigenfunctions are much more symmetric for 9 = 2 (see 
Fig. [3] d and e). This asymmetry has impact on the observed de¬ 
viation between the change of the critical Richardson number for 
9 = 2 and 9 = 10. In fact, in a stratified atmosphere exchanged 
fluid parcels, where one is shifted downwards form the middle 
plane z = 0.5 and the other upwards, move at the same speed 
in opposite directions, but have different Mach numbers. With 
an increasing temperature gradient the effective Mach number 
in the lower half of our domain has a steeper drop such that the 
stabilising effect of compressibility vanishes. Therefore, the sta¬ 
bilising effect of greater Mach numbers in a strongly stratified 
atmosphere is weaker. 


Fig. 4. Both plots show the growth rates £ r in an inviscid weakly ther¬ 
mally stratified atmosphere with 9 = 0.5 for a Richardson number of 
Ri = 0.22. In a) the fluid is compressible with M = 0.09 and in b) 
M = 0.009, which corresponds to a incompressible fluid. The Peclet 
number is varied among three orders of magnitude. 


becomes important in a system where 4 < f B , such that buoyancy 
is much slower than thermal diffusion. As fg has to be smaller 
than the system’s dynamic time scale fg, the Peclet number has 
to become smaller than unity to satisfy these requirements. 

In Fig.[4]the growth rates for different Peclet numbers in a com¬ 
pressible and a weakly compressible fluid are illustrated. Com¬ 
paring instability growth rates for a set up with Pe < 1 and 
Pe > 1 shows that increases with decreasing Pe, where a sig¬ 
nificant jump can be observed when Pe becomes smaller than 
unity. We find that the overall growth rate is smaller for a shear 
flow with Mach number equal to 0.09, this demonstrates a sta¬ 
bilising effect of moderate Mach numbers as discussed in Sect. 
ED- However, an increase in the instability growth rate does 
not necessary indicates a shift of the stability threshold to greater 
Richardson numbers. 


3.2. Small Peclet number regime 

In the following we focus on the impact of thermal diffusion in 
a stably stratified fluid. For non-negligible thermal diffusion the 
non-dimensional Peclet number is given as 


Pe = 


UqL u 
C k ’ 


( 20 ) 


where L„ is the characteristic length of the shear width and Uq 
correspond to the Mach number at the top of the domain as the 
velocity is normalised with respect to the sound speed. The Pe¬ 
clet number is associated with the ratio of advective transport to 
thermal diffusion. 

Varying the Peclet number in an in visci d compressible fluid 
enhances the results found by iLignieres et alJ ( 1999 ) where a 
higher thermal diffusion destabilises the system as it effec¬ 
tively weakens the stable stratification, i.e. the system becomes 
more unstable against buoyancy. Therefore, thermal diffusion 


Focusing on the stability threshold it is necessary to seek for the 
critical Richardson numbers in different Peclet number regimes. 
For previous calculations of the growth rates, the fluid was as¬ 
sumed to be inviscid to simplify the problem and to avoid the 
issue of the initial state being actually not an equilibrium state. 
Here, it is more convenient to include viscosity as it is numeri¬ 
cally easier to obtain results for the limit of small viscosity than 
for an inviscid fluid. 

For a shear flow the onset of instability does not change for 
small enough k, if large enough Reynolds numbers are consid¬ 
ered. However, including viscosity may have non-trivial effects 
on the stability if the fluid has large thermal diffusivity (see l.TonesI 
19771). i ones! (1 977 ) derived a criterion for instability at long 
wavelength of the form kPeRi < 0.086 in an inviscid fluid. For 
a viscous fluid it can be rewritten as kcrReRi < 0.086 such that 
for reasonably small wave numbers, the system can be still un¬ 
stable if the Reynolds number is not sufficiently large. It can be 
partly seen when investigating equation (IT4] > and (fTfil ). where al- 
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most all terms include the wave number, while one of the three 
viscous terms does not. Therefore, for small wave numbers this 
viscous term, which is proportional to the second derivative in 
Z direction, becomes relatively more important for the dynamics 
of the system. To make sure that the results obtained correspond 
to a regime where viscosity does not affect the stability thresh¬ 
old, several computations were repeated with a smaller viscosity 
and a greater viscosity than used in the actual computations. 

In Fig. [5] curves of marginal stability for four different Peclet 
numbers in a polytropic atmosphere are displayed for 9 = 0.5 
and 9 = 2.0. As expected the domain for an unstable shear flow 
increases as the Peclet number is decreased. It shows that not 
only is the growth rate of the instability altered but the stabil¬ 
ity threshold also changes for small Peclet numbers. The critical 
Richardson numbers for Pe < 1 reveal a destabilisation for small 
k. iLignieres et al . ( 1 999) explained this behaviour by the effect 
of the anisotropy of the buoyancy force. The stabilising effect 
of stratification becomes inefficient for predominant horizontal 
motion compared to the thermal diffusion. Indeed, by comput¬ 
ing the ratio of the vertical to horizontal kinetic Energy of the 
unstable mode for a certain k i.e. 

E» = £ w k(z) 2 dz 

Eu £ u k(z) 2 dz 

where w^iz) and Uk(z.) are the eigenfunctions for the vertical and 
horizontal velocity disturbances at a certain k respectively, we 
are able to investigate the nature of the instability. Comparing 
this ratio for a mode with k — 0.1 and k = 0.7 for Pe = 0.1, 
it turns out that the ratio for the larger k-mode is of two orders 
of magnitude greater than for the smaller k-mode. Thus the hor¬ 
izontal motion associated with larger wave lengths in horizontal 
direction is predominant at very small k. 

The results carried out for a greater temperature gradient, 9 = 2 
shown in Fig. [3b, reveal a qualitatively different behaviour for 
the small Peclet numbers than it is the case for 9 = 0.5 shown in 
Fig. 0a. The 9 = 0.5 case has an overall greater critical Richard¬ 
son number, but has a gradual increase towards a smaller critical 
Richardson number in the small k limit. This indicates a more 
efficient destabilisation for most of the wave numbers and a less 
efficient destabilisation for small k. 

Looking again at the ratio of kinetic energy in vertical and hor¬ 
izontal motions reveals that for lower 9 the ratio of vertical mo¬ 
tions to horizontal motions remains significantly smaller than for 
the higher stratification 9 = 2, which indicates that buoyancy 
force is more efficient in a strongly stratified atmosphere against 
the destabilising effect of thermal diffusion. A weakly stratified 
atmosphere can be destabilised quicker by thermal diffusion. 

3.2.1. Effect of the distance to the onset of convection on the 
instability 

The destabilising mechanism of thermal diffusion can be traced 
back to the fact that thermal diffusion weakens the stable strati¬ 
fication against buoyancy. To investigate if, and how, this effect 
changes if the system is further away from the onset of convec¬ 
tion the polytropic index is varied between 1 . 6 , which is close 
to the threshold, and 2.07 which is far from the onset of con¬ 
vection. In Fig. [3 a) the growth rate, A, is plotted for different 
polytropic indices while the Richardson number is fixed to the 
value Ri = 0.22 and the Peclet number is unity. The growth 
rate is given in units of the speed sound, such that a rescal¬ 
ing is necessary while the amplitude of the shear flow is ad¬ 
justed to keep Ri constant. The same is done for a Peclet num- 
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Fig. 5. Plots of critical Richardson numbers for all k in a viscous fluid 
with a Mach number, M = 0.009, a) 9 = 0.5 and b) 9 = 2.0 and viscosity 
(crCk) of order 10 -7 . The Peclet number is varied from 10 to 0.01. 


her of 0.1 in Fig. [6] b) where the same tendency of an increas¬ 
ing growth rate with increasing polytropic index is observed. 
To exclude that this behaviour is due to the effect of thermal 
diffusion the growth rate for two different m is computed for a 
fixed Ri smaller 1 /4 for much greater Pe than unity, where the 
same tendency is found. However, the critical Ri increases as m 
is varied upwards. Therefore, a system with larger polytropic in¬ 
dex, m, is less stable such that for the same Ri the system with 
a larger m is further away from the stability threshold and has 
a higher growth rate. Non-linear direct numerical computations 
of two cases approved the results obtained with the linear EV 
solver. For this purpose equations (Q} - £3 were solved by using 
a hy brid finite-difference/pseudo-spectral cod e (see for exam¬ 
ple iMatthews^eTafUT^ll; IS i I vers et al J l2009al fbL and references 
therein). 

While increasing the polytropic index requires a greater shear 
flow amplitude to obtain the same ratio between the buoyancy 
force and the turnover rate, the available kinetic energy of the 
system is increased. Consequently, as soon as the system can 
overcome the stabilising effect of density stratification the insta¬ 
bility has more available kinetic energy from which it can grow 
more rapidly. 

In the two plots at the bottom of Fig. [3 the growth rates for dif¬ 
ferent polytropic indices are shown while all other parameter are 
held. The solid lines in Fig. [6] c) and d) correspond to the solid 
lines in Fig. [3 a) and b) respectively. As expected the instabil¬ 
ity growth rate decreases with increasing density stratification 
and greater Ri. While for a Peclet number of unity the instability 
shuts down rapidly, which is displayed in Fig. [3c), for Pe < 1 
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Fig. 6. For an inviscid fluid with thermal stratification 6 = 2.0 and a flow of low Mach number of order 10 2 , the growth rate of the linear instability 
is plotted for several parameters. In a) and c) the Peclet number is equal to unity and in b) and d) Pe = 0.1. Ri is varied at the two bottom plots c) 
and d) but it is fixed to Ri = 0.22 for the top plots a) and b). 


the same behaviour is observed, but the instability for greater m 
is still present for small wave numbers. As discussed in Sect. 13.21 
this is explained by the anisotropy of the buoyancy force. 


3.3. Subdominant shear instability 

For certain configurations where the shear width is sufficiently 
small the velocity profile is similar to two counter flows. Then, 
for a stable stratification and Peclet numbers much greater 
than unity a Holmboe like instability is found. The Holm- 
boe instability, which was introduced by tHolmboe dl962h . re¬ 
sults from interacting waves, tha t propagate in opposite direc¬ 
tions dBaines & Mitsuderalll994l) . The Holmboe instability dif¬ 
fers fro m the KH instability in sev eral ways (see the review ar¬ 
ticle bv iPeltier & Caulfield 1200 3). First, the KH instability is 
stationary in a frame of reference, but the Holmboe instability 
has counter propagating unstable modes, which both have the 
same growth rate. One mode occupies the upper plane and the 
other the lower plane, such that a superposition of both solutions 
form a standing wave solution. Second, the Holmboe instabil¬ 
ity is favoured in stable stratified atmospheres and can dominate 
the Kelvin-Helmholtz instability when the stratification stabil¬ 
ity is increased, such thaj_a KH instability disappear due to the 
Richardson criterion (iPeltier & SmvthllT989h . 

In Fig. 0a) there are two arches present. The one with the greater 
maximum growth rate appears at greater k and corresponds to 
the KH instability, where only one unstable mode is present for 
a certain wave number and its phase velocity corresponds to the 
mean flow velocity. The eigenfunction for the vertical velocity of 
these modes is shown in Fig. 0c), which reveals that the insta¬ 
bility is localised at the center of the shear flow. Unstable modes 
with k corresponding to the smaller arch exhibit all properties of 


the Holmboe instability, where two counter propagating modes, 
one in the upper and one in the lower half plane, are present. 
Fig. 0 b) shows the form of the vertical velocity eigenfunction 
that is propagating in the lower half plane, the corresponding 
mode with the opposite phase velocity shows the oscillations in 
the upper half plane. Investigating all possible eigenfunctions at 
the overlap region between the two arches, reveals that Holmboe 
modes are present but have a smaller growth rate as the KH in¬ 
stability or vice versa. 

The observed Holmboe like instability only appears for large Pe¬ 
clet numbers while for small Pe the range for the KH instabil¬ 
ity is enlarged to smaller wave numbers, such that even in the 
presents of the Holmboe like instability the KH instability for the 
same wave numbers dominates. This can be seen in Fig. 0a) and 
b), where a small shoulder is present in the range 0.35 < k < 0.5 
in a) for Pe — 1 is dominated by the KH instability in b) for 
Pe = 0.1. 

To make sure that these modes are physical and not a numeri¬ 
cal artifact we checked that they remain for a system where we 
assume an inviscid non diffusive incompressible fluid. In a sec¬ 
ond check, we consider a non-uniform grid distribution in the z 
direction, with more resolution where the largest gradients are 
observed. With both approaches a second instability, with the 
same properties, was found at slightly smaller wave numbers, 
k, than the KH instability. Using non-linear computations both 
instabilities are found during the linear growth regime with sim¬ 
ilar growth rates and eigenfunctions, as predicted from the linear 
computations. 

4. Conclusions 

Shear flow instabilities, that may lead to turbulent motions, are 
important for the understanding of dynamics in stellar interiors 
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Fig. 7. The growth rate for an inviscid fluid with a temperature gradient, 
6 = 1.0, Peclet number of order 10 3 and Ri = 0.3 is shown in plot a). 
The form of the eigenfunctions found for the vertical velocity w in the 
small arch and the large arch are displayed in b) and c) respectively. 


and in other astrophysical flows. Using linear stability analysis 
the onset of shear flow instability in a compressible polytropic 
atmosphere was investigated, where a detailed analysis of the ef¬ 
fect of moderate Mach numbers, small Peclet numbers and dif¬ 
ferent polytropic indices was carried out separately. 

For a flow of moderate Mach numbers we found a stabilising 
effect that results in a critical Richardson number which is less 
than 1/4, such that for a certain range of Ri less than 1/4 the 
shear flow remains stable to KH instability. This effect becomes 
significant for Mach numbers greater than 0.1. However, the sta¬ 
bilising is weaker for strongly stratified atmospheres. As this 
effect is most relevant at moderate to high Mach numbers it 
can be neglected for most stellar regions, in particular for the 
tachocline. However, this result might be important in other as¬ 
trophysical objects where high Mach numbers are common. 

For fluids with high thermal diffusivity, where the Peclet number 
drops below unity, a destabilisation of the system can be shown. 
Our results are in agreement with lLignieres et aTI (11999b . where 
the regime of small Peclet numbers were examined in a Boussi- 
nesq fluid. We find a significantly greater growth rate for highly 
diffusive fluids as well as greater critical Richardson numbers for 
Pe < 1 compared to systems with large Peclet numbers. Gener¬ 
ally, steeper temperature gradients lead to an overall stabilisation 
for small Peclet numbers. However, it should be noted that, for 
very small wave numbers we found that the opposite is true due 
to anisotropy of the buoyancy force. 

An interesting result is that we find the possibility of Holmboe 
like instabilities present in polytropic atmopheres. This instabil¬ 
ity is dominated by the KH instability for small Peclet numbers, 
but is clearly present for large Peclet numbers. While no setup 
was found where the Holmboe like instability dominates the KH 
instability it may be found in future investigations. 

Studying more complex systems, which include key properties 
of stars, by means of a linear stability analysis provides a power¬ 
ful tool to investigate their linear behaviour. However, the verti¬ 
cal length scale for which these shear flows can become unstable 
is far below the resolution of helioseismological techniques and 


as th e perturbation’s ampli tude can be finite in the tachocline 
(see lMichaud & Zahnll998l) the assumption of infinitesimal per¬ 
turbations made in the stability analysis might lead to different 
results. To clarify the second issue it is necessary to investigate 
the stability properties of such flows by means of direct numer¬ 
ical computations, where the initial state is perturbed by pertur¬ 
bations with a finite amplitude. In addition to finite amplitude 
effects, the non-linear behaviour of unstable shear flows after a 
saturation is crucially important for the understanding of mixing 
in stellar interiors regardless of the nature of the initial instabil¬ 
ity. 

Having established an understanding of the hydrodynamical 
problem, further non-linear investigations are underway to ob¬ 
tain a full picture. 
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